
/*************
author:: libo
license:: MIT
*****************/

#include "Mesh/Cell_Traits.h"
#include "Mesh/_Func_.h"
#include "clapack/include/f2c.h"
#include "tools/intrbtree.h"
#include "tools_node.h"
#include<tool/libcell_tools_algorithm.h>

#include<Mesh/libcell_iterator.h>
#include <Math/LB_Math.h>
#include <Math/LB_Quaternions.h>
#include <Matrix/LB_Matrix.h>
#include<Algorithm/threeD_elementary_geometry_tool_library.h>
//#include <unordered_map>
#include<tools/inthashmap.h>
#include<tools_vector.h>

#define quote lib_cell_quote


double* mesh_compute_bounding_box_from_node(Node* n)
{
    if(n==NULL)
    {
        return NULL;
    }
    double* re=(double*)malloc(sizeof(double)*6);
    template_v* v=(template_v*)(n->value);
    re[0]=v->point[0];re[1]=v->point[1];re[2]=v->point[2];
    re[3]=re[0];re[4]=re[1];re[5]=re[2];
    for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
    {
        v=(template_v*)(nit->value);
        for(int i=0;i<3;i++)
        {
            v->point[i]>re[3+i]?(re[3+i]=v->point[i]):(re[i+3]=re[i+3]);
            v->point[i]<re[i]?(re[i]=v->point[i]):(re[i]=re[i]);
        }
    }
    return re;
}
double* mesh_compute_bounding_box(Mesh* m)
{
    double * re=(double*)malloc(sizeof(double)*6);
    auto it=m->v_begin(m);
    template_v* v=quote(it);
    re[0]=v->point[0];re[1]=v->point[1];re[2]=v->point[2];
    re[3]=re[0];re[4]=re[1];re[5]=re[2];
    it++;
    for(;it!=m->v_end(m);it++)
    {
        v=quote(it);
        for(int i=0;i<3;i++)
        {
            v->point[i]>re[3+i]?(re[3+i]=v->point[i]):(re[i+3]=re[i+3]);
            v->point[i]<re[i]?(re[i]=v->point[i]):(re[i]=re[i]);
        }
    }
    return re;
}
void mesh_adjust_mesh_to_unit_hexahedron(Mesh* m,double* box)
{
    if(box==NULL){
        // double *box1=mesh_compute_bounding_box(m);
        // if(box1==NULL){return ;}
        // memmove(box,box1,sizeof(double)*6);
        // SAFE_FREE(box1);
        return ;
    }

    //printf("box:%lf %lf %lf %lf %lf %lf\n",box[0],box[1],box[2],box[3],box[4],box[5]);
    double centroid[3]={(box[0]+box[3])/2.0,(box[1]+box[4])/2.0,(box[2]+box[5])/2.0};
    //printf("centroid : %lf %lf %lf\n",centroid[0],centroid[1],centroid[2]);
    double max=box[3]-box[0];

    max=((box[4]-box[1])>max?((box[4]-box[1])):(max));
    max=((box[5]-box[2])>max?((box[5]-box[2])):(max));
    // printf("max:%lf\n",max);
    for(auto vit=m->v_begin(m);vit!=m->v_end(m);vit++)
    {
        template_v*v=(template_v*)(quote(vit));
        v->point[0]-=centroid[0];
        v->point[1]-=centroid[1];
        v->point[2]-=centroid[2];
        v->point[0]*= (2.0/max);
        v->point[1]*= (2.0/max);
        v->point[2]*= (2.0/max);

    }
}

// static int is_two_array_equal_ignoring_order(template_v** vs1,template_v** vs2,int len)
// {
//     for(int i=0;i<len;i++)
//     {
//         int flag=0;
//         for(int j=0;j<len;j++)
//         {
//             if(vs1[i]==vs2[j]){flag=1;break;}
//         }
//         if(flag==0){return 0;}
//     }
//     return 1;
// }


Node* Mesh_adjust_halffaces(template_c* c)
{

    Node* hfs=NULL;
    for(int i=c->halffaces_size-1;i>=0;i--)
    {
        hfs=node_overlying(hfs,c->halffaces[i]);
    }
    //node_copy(c->halffaces);

    Node* re=NULL;
    re=node_pushback(re,hfs->value);
    hfs=node_remove(hfs,hfs);
    int flag=1;
    while(flag)
    {
        template_hf* hf=(template_hf*)(re->value);
        if(hf->vertices_size<=1)
        {
            free_node(hfs);
            return re;
        }
        template_v* v=hf->vertices[1];
        flag=0;
        for(Node* nit=hfs;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_hf* hf1=(template_hf*)(nit->value);
            if(Mesh_vertex_is_in_face(v,hf1->face))
            {
                flag=1;
                re=node_pushback(re,nit->value);
                hfs=node_remove(hfs,nit);
                break;
            }
        }
    }
    free_node(hfs);

    return node_reverse(re);

}
template_v** mesh_adjust_vertices(Mesh*m,template_c*c)
{
    Node*node= Mesh_adjust_halffaces(c);
    template_v**re=(template_v**)malloc(sizeof(template_v*)*c->vertices_size);
    int i=0;
    for(Node*nit=node;nit->Next!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf*hf=(template_hf*)(nit->value);
        re[i]=hf->vertices[1];
        i++;
    }
    free_node(node);

    return re;
}


static Node* split_unicom_region_from_sub_node_with_mark_faces(LB_Vector1_voidptr* vec ,Int_RB_Tree*tree_f )
{
    Node* re=NULL;

    Int_RB_Tree* tree0=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree0);

    // for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    for(int i=0;i<vec->size;i++)
    {
        template_c*c=(template_c*)(vec->extra_data[i]);
        // template_c* c=(template_c*)(nit->value);
        tree0->insert(tree0,c->id,c);
    }

    while(tree0->size>0)
    {
        template_c* c=(template_c* )((tree0->begin(tree0)).second);
        Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
        int_rb_tree_init(tree);
        tree->insert(tree,c->id,c);
        tree0->erase(tree0,c->id);
        //*****************
        Node* temp=NULL;
        temp=node_overlying(temp,c);

        Node* temp1=NULL;
        while(temp!=NULL)
        {
            for(Node* nit=temp;nit!=NULL;nit=(Node*)(nit->Next))
            {
                c=(template_c*)(nit->value);
                //for(Node*nit1=c->halffaces;nit1!=NULL;nit1=(Node*)(nit1->Next) )
                for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++ )
                {
                    if(tree_f->find(tree_f, quote(chfit)->face->id)!=NULL)
                    {
                        continue;
                    }
                    template_c* c3=Mesh_s_opposite_halfface(quote(chfit))->cell;

                    if(c3==NULL||tree->find(tree,c3->id)!=NULL||tree0->find(tree0,c3->id)==NULL )
                    {
                        continue;
                    }
                    else
                    {
                        tree->insert(tree,c3->id,c3);
                        tree0->erase(tree0,c3->id);
                        temp1=node_overlying(temp1,c3);
                    }
                }
            }
            free_node(temp);
            temp=temp1;
            temp1=NULL;
        }
        //printf("%d\n",tree->size);
        //******************
        re=node_overlying(re,tree);
    }
    int_rb_tree_free(tree0);
    return re;
}

static template_c* my_get_cells_copy_data(template_c* c)
{
    template_c*re=(template_c*)malloc(sizeof(template_c));
    Cell_init_(re);
    re->id=c->id;
    //for(Node* nit=c->halffaces;nit!=NULL;nit=(Node*)(nit->Next))
    re->halffaces_size=c->halffaces_size;
    re->halffaces=(template_hf**)malloc(sizeof(template_hf*)*re->halffaces_size);
    int i=0;
    for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++)
    {
        template_hf* hf=quote(chfit);
        template_hf*hf1 =(template_hf*)malloc(sizeof(template_hf));
        HalfFace_init_(hf1);
        //re->halffaces=node_pushback(re->halffaces,hf1);
        re->halffaces[i]=hf1;
        hf1->vertices_size=hf->vertices_size;
        hf1->vertices=(template_v**)malloc(sizeof(template_v*)*hf->vertices_size);
        memmove(hf1->vertices,hf->vertices,sizeof(template_v*)*hf->vertices_size);
        i++;
    }

    return re;
}
static void my_free_copy_cell_data(template_c*cc)
{
    for(auto chfit=Mesh_chf_begin(NULL,*cc);chfit!=Mesh_chf_end(NULL,*cc);chfit++)
    {
        template_hf*hf=quote(chfit);
        free_HalfFace(hf);
        // SAFE_FREE(hf->vertices);
        // free(hf);
    }

    free_Cell(cc);
}

static template_c* my_create_copy_data_cell(Mesh* m,template_c*cc)
{
    //int len=node_size(cc->halffaces);
    int len=cc->halffaces_size;
    template_hf**hfs=(template_hf**)malloc(sizeof(template_hf*)*len);
    //int i=0;
    for(int i=0;i<len;i++)
    {
        template_hf* hf=cc->halffaces[i];
        template_f*f=m->create_facev(m,hf->vertices,hf->vertices_size);
        hfs[i]=m->create_halfface(m,f,hf->vertices,hf->vertices_size);
    }
    template_c*re=m->create_cellf(m,hfs,len);
    free(hfs);
    return re;
}

static Node*  my_tool_copy_new_cells(Mesh* m,template_v*v,Int_RB_Tree*cells ,Int_RB_Tree*re,Int_RB_Tree *re1)
{
    Node* re_node=NULL;
    template_v* nv=m->create_vertexv(m,v->point,3);
    //printf("cuda0\n");
    re1->insert(re1,nv->id,v);
    Node* temp_node=NULL;
    for(auto it=cells->begin(cells);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        template_c*cc= my_get_cells_copy_data(c);
        //for(Node* nit=cc->halffaces;nit!=NULL;nit=(Node*)(nit->Next))
        for(auto chfit=m->chf_begin(m,*cc);chfit!=m->chf_end(m,*cc);chfit++)
        {
            template_hf*hf=quote(chfit);
            for(int j=0;j<hf->vertices_size;j++)
            {
                if(hf->vertices[j]==v)
                {
                    hf->vertices[j]=nv;
                    break;
                }
            }
        }
        temp_node=node_overlying(temp_node,cc);
    }
    //printf("cuda2\n");
    for(auto it=cells->begin(cells);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        m->delete_cell(m,*c,true);
        // template_c*c=m->get_cellp(m,it.first);
        // if(c!=NULL)
        // {
        //     m->delete_cell(m,*c,true);
        // }
    }
    //printf("cuda1\n");
    for(Node* nit=temp_node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c*cc=(template_c*)(nit->value);

        template_c*c=my_create_copy_data_cell(m,cc);
        if(c!=NULL)
        {
            re_node=node_overlying(re_node,c);
            if(re!=NULL)
            {
                int*id_value=(int*)malloc(sizeof(int));
                *id_value=cc->id;
                re->insert(re,c->id,id_value);
            }
        }


        my_free_copy_cell_data(cc);
    }

    free_node(temp_node);
    return re_node;
}

void Mesh_cut_along_the_curvef(Mesh* m,Node* faces,Int_RB_Tree* re,Int_RB_Tree* re1)
{
    //printf("mesh cut along_the curve f node :%d\n",node_size(faces) );
    m->manifold_require=1;
    Int_RB_Tree* tree_f=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_f);
    Int_RB_Tree* tree_fv=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_fv);
    for(Node* nit=faces;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_f*f=(template_f*)(nit->value);
        if(f==NULL){continue;}
        tree_f->insert(tree_f,f->id,f);
        tree_fv->insert(tree_fv,f->vertices[0]->id,f->vertices[0]);
        tree_fv->insert(tree_fv,f->vertices[1]->id,f->vertices[1]);
    }
    printf("here1 %d\n",tree_f->size);
    for(auto it=tree_fv->begin(tree_fv);it.it!=NULL;it++)
    {
        template_v* v=(template_v*)(it.second);
        // printf("once1  ");
        Node* split_cells=split_unicom_region_from_sub_node_with_mark_faces(&(v->cellsn),tree_f);

        if(split_cells==NULL){continue;}
        // printf("once4\n");
        for(Node* nit=(Node*)(split_cells->Next);nit!=NULL;nit=(Node*)(nit->Next))
        {
            Int_RB_Tree* tree=(Int_RB_Tree*)(nit->value);
            Node* temp=my_tool_copy_new_cells(m,v,tree,re,re1);

            free_node(temp);
            int_rb_tree_free(tree);
        }
        //printf("once2\n");
        int_rb_tree_free((Int_RB_Tree*)(split_cells->value));
        free_node(split_cells);
        //printf("once3\n");
    }

    int_rb_tree_free(tree_f);
    int_rb_tree_free(tree_fv);

    return ;
}

Node* split_unicom_region_from_sub_node( Node*node )
{
    Node* re=NULL;
    if(node==NULL){return re;}
    Int_RB_Tree* tree0=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree0);
    for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        tree0->insert(tree0,c->id,c);
    }

    while(tree0->size>0)
    {
        template_c* c=(template_c* )(tree0->begin(tree0).second);
        Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
        int_rb_tree_init(tree);
        tree->insert(tree,c->id,c);
        tree0->erase(tree0,c->id);
        //*****************
        Node* temp=NULL;
        temp=node_overlying(temp,c);

        Node* temp1=NULL;
        while(temp!=NULL)
        {
            for(Node* nit=temp;nit!=NULL;nit=(Node*)(nit->Next))
            {
                c=(template_c*)(nit->value);
                //for(Node*nit1=c->halffaces;nit1!=NULL;nit1=(Node*)(nit1->Next) )
                for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++)
                {
                    template_c* c3=Mesh_s_opposite_halfface(quote(chfit))->cell;
                    if(c3==NULL||tree->find(tree,c3->id)!=NULL||tree0->find(tree0,c3->id)==NULL )
                    {
                        continue;
                    }
                    else
                    {
                        tree->insert(tree,c3->id,c3);
                        tree0->erase(tree0,c3->id);
                        temp1=node_overlying(temp1,c3);
                    }
                }
            }
            free_node(temp);
            temp=temp1;
            temp1=NULL;
        }
        //printf("%d\n",tree->size);
        //******************
        re=node_overlying(re,tree);
    }

    int_rb_tree_free(tree0);
    return re;
}

Node* split_unicom_region_from_sub_vector( LB_Vector1_voidptr vec )
{
    Node* re=NULL;
    if(vec.size<=0){return re;}
    Int_RB_Tree* tree0=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree0);
    // for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    for(int i=0;i<vec.size;i++)
    {
        // template_c* c=(template_c*)(nit->value);
        template_c*c=(template_c*)(*lb_vector1_at_voidptr(
            &(vec), i));

        tree0->insert(tree0,c->id,c);
    }

    while(tree0->size>0)
    {
        template_c* c=(template_c* )(tree0->begin(tree0).second);
        Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
        int_rb_tree_init(tree);
        tree->insert(tree,c->id,c);
        tree0->erase(tree0,c->id);
        //*****************
        Node* temp=NULL;
        temp=node_overlying(temp,c);

        Node* temp1=NULL;
        while(temp!=NULL)
        {
            for(Node* nit=temp;nit!=NULL;nit=(Node*)(nit->Next))
            {
                c=(template_c*)(nit->value);
                //for(Node*nit1=c->halffaces;nit1!=NULL;nit1=(Node*)(nit1->Next) )
                for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++)
                {
                    template_c* c3=Mesh_s_opposite_halfface(quote(chfit))->cell;
                    if(c3==NULL||tree->find(tree,c3->id)!=NULL||tree0->find(tree0,c3->id)==NULL )
                    {
                        continue;
                    }
                    else
                    {
                        tree->insert(tree,c3->id,c3);
                        tree0->erase(tree0,c3->id);
                        temp1=node_overlying(temp1,c3);
                    }
                }
            }
            free_node(temp);
            temp=temp1;
            temp1=NULL;
        }
        //printf("%d\n",tree->size);
        //******************
        re=node_overlying(re,tree);
    }

    int_rb_tree_free(tree0);
    return re;

}

Node* split_unicom_region( Int_RB_Tree* tree_c )
{
    Node* re=NULL;
    if(tree_c==NULL||tree_c->size<=0){return re;}
    Int_RB_Tree* tree0=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree0);
    // Int_Hash_Map *hmap=(Int_Hash_Map*)malloc(sizeof(Int_Hash_Map));
    // int_hash_map_init(hmap,50);
    // std::unordered_map<int, template_c*> tree0;
    void** temp_pro=(void**)malloc(sizeof(void*)*tree_c->size);

    int i=0;
    int mark0=0,mark1=1;
    for(auto it= tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        temp_pro[i]=c->prop;
        c->prop=&mark0;
        // hmap->insert(hmap,c->id,c);
        tree0->insert(tree0,c->id,c);
        //tree0.emplace(c->id, c);
        i++;
    }


    LB_List ll1,ll2;
    lb_list_init(&ll1);lb_list_init(&ll2);
    LB_List  *llp1=&ll1,*llp2=&ll2;

    while(tree0->size>0)
    {
        //template_c* c=(template_c* )(hmap->begin(hmap).second);
        template_c* c=(template_c* )(tree0->begin(tree0).second);
        // auto it = tree0.begin();
        // template_c* c = it->second;
        Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
        int_rb_tree_init(tree);
        tree->insert(tree,c->id,c);
        //hmap->erase(hmap,c->id);
        tree0->erase(tree0,c->id);
        //tree0.erase(it);
        c->prop=&mark1;
        //*****************
        // Node* temp=NULL;
        // temp=node_overlying(temp,c);
        // vecp->size=0;
        // vecp->push_back(vecp,c);
        llp1->size=0;llp1->end=NULL;
        lb_list_push_back(llp1,c);
        //Node* temp1=NULL;
        llp2->size=0;llp2->end=NULL;
        //vecp1->size=0;
        //
        //while(vecp->size>0)
        //while(temp!=NULL)
        while(llp1->end!=NULL)
        {
            //
            //for(int i=0;i<vecp->size;i++)
            //for(Node* nit=temp;nit!=NULL;nit=(Node*)(nit->Next))
            for(Node* nit=llp1->end;nit!=NULL;nit=(Node*)(nit->Prev))
            {
                //c=(template_c*)((vecp->at(vecp,i))[0]);
                c=(template_c*)(nit->value);
                //for(Node* chfit=c->halffaces; chfit!=NULL;chfit=(Node*)(chfit->Next))
                for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++)
                {
                    template_c* c3=Mesh_s_opposite_halfface(quote(chfit))->cell;
                    if(c3==NULL||c3->prop==&mark1||( c3->prop!=&mark1&&c3->prop!=&mark0 ) )
                    {
                        continue;
                    }
                    else
                    {
                        tree->insert(tree,c3->id,c3);
                        c3->prop=&mark1;
                        //hmap->erase(hmap,c3->id);
                        tree0->erase(tree0,c3->id);
                        //tree0.erase(c3->id);
                        //temp1=node_overlying(temp1,c3);
                        lb_list_push_back(llp2,c3);
                        //vecp1->push_back(vecp1,c3);
                    }
                }
            }
            // LB_Vector_voidptr* tempvec=vecp;
            // vecp=vecp1;
            // vecp1=tempvec;
            // vecp1->size=0;
            // free_node(temp);
            // temp=temp1;
            // temp1=NULL;
            llp1->size=0;llp1->end=NULL;
            LB_List * templl=llp1;
            llp1=llp2;
            llp2=templl;
        }
        printf("%d\n",tree->size);
        //******************
        re=node_overlying(re,tree);
    }
    i=0;
    for(auto it= tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        c->prop= temp_pro[i];
        i++;
    }
    free(temp_pro);
    //int_hash_map_free(hmap);
    int_rb_tree_free(tree0);
    lb_list_clear(&ll1);lb_list_clear(&ll2);
    //vec.clear(&vec);vec1.clear(&vec1);
    return re;
}

Node* split_unicom_region2mesh(Mesh *m )
{
    Node* re=NULL;
    Node* node=split_unicom_region(m->cells);
    if(node==NULL){return re;}
    for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        Int_RB_Tree* tree=(Int_RB_Tree*)(nit->value);
        Mesh * mesh=(Mesh*)malloc(sizeof(Mesh));
        Mesh_init(mesh);
        mesh->dimension=m->dimension;
        mesh->simplex=m->simplex;
        mesh->manifold_require=m->manifold_require;
//*********************

        Int_RB_Tree* treev=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
        int_rb_tree_init(treev);

        for(auto it=tree->begin(tree);it.it!=NULL;it++)
        {
            template_c* c=(template_c*)(it.second);
            for(int j=0;j<c->vertices_size;j++)
            {
                if(treev->find(treev,c->vertices[j]->id)==NULL)
                {
                    template_v*v= mesh->create_vertexv(mesh,c->vertices[j]->point,c->vertices[j]->point_size);

                    treev->insert(treev,c->vertices[j]->id,v);
                }
            }

            int j=0;
            //int hfs_l=node_size(c->halffaces);
            int hfs_l=c->halffaces_size;
            template_hf**hfs=(template_hf**)malloc(sizeof(template_hf*)*hfs_l);
            for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
            {
                template_v**vs=(template_v**)malloc(sizeof(template_v*)*quote(chfit)->vertices_size);
                for(int i=0;i<quote(chfit)->vertices_size;i++)
                {
                    vs[i]=(template_v*)(treev->find(treev,quote(chfit)->vertices[i]->id));
                }

                template_f*f=mesh->create_facev(mesh, vs,quote(chfit)->vertices_size );
                hfs[j]=mesh->create_halfface(mesh,f,vs,quote(chfit)->vertices_size);
                j++;
                free(vs);
            }

            mesh->create_cellf(mesh,hfs,hfs_l);
            free(hfs);
        }

        int_rb_tree_free(treev);
//*******************
        re=node_overlying(re,mesh);

        int_rb_tree_free(tree);
    }

    free_node(node);

    return re;
}

static template_v* get_next_shortest_vertex(Mesh*m,Int_RB_Tree* tree,Int_RB_Tree* tree1,Node* vs)
{
    template_v* v=(template_v*)(vs->value);
    //if(v==NULL){printf("cuowwwww\n");}
    for(auto vfit=m->vf_begin(m,*v);vfit!=m->vf_end(m,*v);vfit++)
    {
        template_v* v1=NULL;
        (quote(vfit)->vertices[0]==v)?(v1=quote(vfit)->vertices[1]):(v1=quote(vfit)->vertices[0]);
        //if(v==v1){printf("cuow\n");}
        double* value=(double*)(tree->find(tree,v1->id));
        double* value1=(double*)(tree->find(tree,v->id));
        //if(value1==NULL){printf("cuowu\n");}
        if(value==NULL)
        {
            value=(double*)malloc(sizeof(double));
            *value=*value1+ distance_of_two_points(v->point,v1->point,3);
            tree->insert(tree,v1->id,value);
            //printf("inset:%d %lf\n",v1->id,*value);
            tree1->insert(tree1,v1->id,value);
        }
        else
        {
            if(*value>*value1+distance_of_two_points(v->point,v1->point,3))
            {
                *value=*value1+ distance_of_two_points(v->point,v1->point,3);
            }
        }
        //printf("once\n");

    }
    //printf("here\n");

    //printf("tree1->size:%d\n",tree1->size);
    auto it=tree1->begin(tree1);
    int id=it.first;
    if(it.second==NULL){return NULL;}
    double temp_value=*((double*)(it.second));

    for(;it.it!=NULL;it++)
    {
        if(*((double*)(it.second))<temp_value)
        {
            id=it.first;
            temp_value=*((double*)(it.second));
        }
    }
    //printf("here1 %lf %d\n",temp_value,id);

    return m->get_vertexp(m,id);
}
static Node* get_shortest_path1(Mesh* m, template_v* v1,template_v*v2,Int_RB_Tree* tree)
{
    Node* re=NULL;
    template_v* v=v2;
    while(v!=v1)
    {
        re=node_overlying(re,v);
        //printf("ud:%d\n",v->id);
        template_v* v3=NULL,*v4=NULL;
        double temp=0;

        for(auto vfit=m->vf_begin(m,*v);vfit!=m->vf_end(m,*v);vfit++)
        {
            (quote(vfit)->vertices[0]==v)?(v3=quote(vfit)->vertices[1]):(v3=quote(vfit)->vertices[0]);
            double *value=(double*)(tree->find(tree,v3->id));
            if(value==NULL)
            {
                continue;
            }

            if(*value<temp||v4==NULL)
            {
                temp=*value;
                v4=v3;
            }
        }
        v=v4;
    }
    re=node_overlying(re,v1);

    return re;
}

Node* get_shortest_path(Mesh* m, template_v* v1,template_v*v2)
{
    if(v1==NULL||v2==NULL){return NULL;}
    Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree);

    Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    Node* vs=NULL;
    vs=node_overlying(vs,v1);

    double*value=(double*)malloc(sizeof(double));
    *value=0;
    tree->insert(tree,v1->id,value);

    //printf("begin shortest \n");
    template_v *v=NULL;
    //int nums=0;
    do
    {
        v =get_next_shortest_vertex(m,tree,tree1,vs);
        if(v==NULL){
            printf("isolate vertex\n");
         break;}
        vs=node_overlying(vs,v);
        tree1->erase(tree1,v->id);
        //printf("erase:%d %d\n",v->id,v2->id);
        // if(nums>60)
        // {
        //     break;
        // }
        // nums++;
    }while(v!=v2);
    // printf("begin shortest \n");

    Node* re1=NULL;
    if(v!=NULL)
    {
        re1=get_shortest_path1(m, v1,v2,tree);
    }


    //printf("end shortest \n");

    int_rb_tree_free_value(tree);free_node(vs);
    int_rb_tree_free(tree);int_rb_tree_free(tree1);

    return re1;
}

double compute_surface_area(Mesh *m)
{
    double re=0;
    for(auto cit=m->c_begin(m);cit!=m->c_end(m);cit++)
    {
        template_v* v=quote(cit)->vertices[0];
        for(auto chfit=m->chf_begin(m,*cit);chfit!=m->chf_end(m,*cit);chfit++)
        {
            if(quote(chfit)->vertices[0]!=v&&quote(chfit)->vertices[1]!=v )
            {

                re+=compute_area_of_triangle(v->point,quote(chfit)->vertices[0]->point,
                    quote(chfit)->vertices[1]->point );
            }
        }
    }

    return re;
}

static inline void  get_cube_from_cell(double **ps,int len,double *cube)
{
    //memset(cube,0,sizeof(double)*4);
    cube[0] = ps[0][0];    cube[1] = ps[0][1];
    cube[2] = ps[0][0];    cube[3] = ps[0][1];
    for(int i=1;i<len;i++)
    {
        ps[i][0]<cube[0]?(cube[0]=ps[i][0]):(ps[i][0]>cube[2]?(cube[2]=ps[i][0]):(cube[2]=cube[2]));
        ps[i][1]<cube[1]?(cube[1]=ps[i][1]):(ps[i][1]>cube[3]?(cube[3]=ps[i][1]):(cube[3]=cube[3]));

    }
    // p2[0]<cube[0]?(cube[0]=p2[0]):(p2[0]>cube[2]?(cube[2]=p2[0]):(cube[2]=cube[2]));
    // p2[1]<cube[1]?(cube[1]=p2[1]):(p2[1]>cube[3]?(cube[3]=p2[1]):(cube[3]=cube[3]));

    // p3[0]<cube[0]?(cube[0]=p3[0]):(p3[0]>cube[2]?(cube[2]=p3[0]):(cube[2]=cube[2]));
    // p3[1]<cube[1]?(cube[1]=p3[1]):(p3[1]>cube[3]?(cube[3]=p3[1]):(cube[3]=cube[3]));

}
static inline void fill_center_data(template_c* c,double* cen)
{
    double min[3]={c->vertices[0]->point[0],c->vertices[0]->point[1],c->vertices[0]->point[2]};
    double max[3]={min[0],min[1],min[2]};
    for(int i=1;i<c->vertices_size;i++)
    {
        c->vertices[i]->point[0]<min[0]?(min[0]=c->vertices[i]->point[0]):(
                c->vertices[i]->point[0]>max[0]?(max[0]=c->vertices[i]->point[0]):(max[0]=max[0]));

        c->vertices[i]->point[1]<min[1]?(min[1]=c->vertices[i]->point[1]):(
                c->vertices[i]->point[1]>max[1]?(max[1]=c->vertices[i]->point[1]):(max[1]=max[1]));

        c->vertices[i]->point[2]<min[2]?(min[2]=c->vertices[i]->point[2]):(
                c->vertices[i]->point[2]>max[2]?(max[2]=c->vertices[i]->point[2]):(max[2]=max[2]));

    }
    cen[0]=(min[0]+max[0])/2.0;cen[1]=(min[1]+max[1])/2.0;cen[2]=(min[2]+max[2])/2.0;
}

static void sampling_from_one_cell(Mesh* m,template_c*c,Node**n,double dense )
{
    //printf("cid :%d %d\n",c->id,c->vertices_size);
    double cen[3]={0};
    fill_center_data(c,cen);

    double** ps=(double**)malloc(sizeof(double*)*c->vertices_size);
    for(int i=0;i<c->vertices_size;i++)
    {
        ps[i]=(double*)malloc(sizeof(double)*3);
        ps[i][0]=c->vertices[i]->point[0]-cen[0];
        ps[i][1]=c->vertices[i]->point[1]-cen[1];
        ps[i][2]=c->vertices[i]->point[2]-cen[2];
    }

    double q1[3]={ps[1][0]-ps[0][0],ps[1][1]-ps[0][1],ps[1][2]-ps[0][2]};
    double q2[3]={ps[2][0]-ps[0][0],ps[2][1]-ps[0][1],ps[2][2]-ps[0][2]};

    double *normal=out_product(q1,q2);
    if(!normalize(normal,3))
    {
        SAFE_FREE(normal);
        return ;
    }
    double theta = acos(normal[2]);
    double  n1[3]={normal[1],-normal[0],0};
    if(normalize(n1,3))
    {
        double*matd = get_matrix_from_rotation(theta,n1);

        for(int k=0;k<c->vertices_size;k++)
        {
            double pp[3]={0};
            for(int i=0;i<3;i++)
            {
                for(int j=0;j<3;j++)
                {
                    pp[i]+=(matd[i*3+j]*ps[k][j]);

                }
            }
            //printf("pp:%lf %lf %lf\n",pp[0],pp[1],pp[2]);
            memmove(ps[k],pp,sizeof(double)*3);
        }
        free(matd);
    }

    double cube[4]={0};
    get_cube_from_cell(ps,c->vertices_size,cube);

    int min_x=(int)(floor(cube[0]/dense));
    int min_y=(int)(floor(cube[1]/dense));

    int max_x=(int)(ceil(cube[2]/dense));
    int max_y=(int)(ceil(cube[3]/dense));
    LB_Matrix** mats=(LB_Matrix**)malloc(sizeof(LB_Matrix*)*(c->vertices_size-2));
    LB_Matrix** matis=(LB_Matrix**)malloc(sizeof(LB_Matrix*)*(c->vertices_size-2));
    // printf("here1\n");
    for(int i=0;i<c->vertices_size-2;i++)
    {
        mats[i]=(LB_Matrix*)malloc(sizeof(LB_Matrix));
        lb_matrix_init_double(mats[i],2,2);
        double* data=(double*)(mats[i]->data);
        data[0]=ps[i+1][0]-ps[0][0]; data[1]=ps[i+2][0]-ps[0][0];
        data[2]=ps[i+1][1]-ps[0][1]; data[3]=ps[i+2][1]-ps[0][1];
        matis[i]=mats[i]->inverse(mats[i]);

    }
    //printf("here2\n");

    // LB_Matrix* mati=mat->inverse(mat);
    // data=(double*)(mati->data);
    double y[3]={0},pp[2]={0};
    for(int i=min_x;i<max_x;i++)
    {
        pp[0]=(i+(rand()%30000)/30000.0)*dense-ps[0][0];

        for(int j=min_y;j<max_y;j++)
        {
            pp[1]=(j+(rand()%30000)/30000.0)*dense-ps[0][1];
            for(int k=0;k<c->vertices_size-2;k++)
            {
                if(matis[k]==NULL){continue;}
                double *data=(double*)(matis[k]->data);
                y[1]=pp[0]*data[0]+pp[1]*data[1];
                y[2]=pp[0]*data[2]+pp[1]*data[3];
                y[0]=1-y[1]-y[2];
                if(y[0]>0&&y[1]>0&&y[2]>0)
                {
                    Vertex*v=(Vertex*)malloc(sizeof(Vertex));
                    Vertex_init_(v);
                    v->point_size=3;
                    v->point=(double*)malloc(sizeof(double)*3);
                    v->point[0]=y[0]*c->vertices[0]->point[0]+y[1]*c->vertices[k+1]->point[0]+y[2]*c->vertices[k+2]->point[0];
                    v->point[1]=y[0]*c->vertices[0]->point[1]+y[1]*c->vertices[k+1]->point[1]+y[2]*c->vertices[k+2]->point[1];
                    v->point[2]=y[0]*c->vertices[0]->point[2]+y[1]*c->vertices[k+1]->point[2]+y[2]*c->vertices[k+2]->point[2];
                    v->prop=c;
                    *n =node_overlying(*n,v);
                    break;
                }
            }
        }
    }

    for(int i=0;i<c->vertices_size-2;i++)
    {
        lb_matrix_free(mats[i]);
        lb_matrix_free(matis[i]);
    }

    free(mats);free(matis);
    for(int i=0;i<c->vertices_size;i++)
    {
        free(ps[i]);
    }
    free(ps);
}
Node* create_surface_uniform_sampling_points(Mesh* m,int num,double dense)
{
    Node* re=NULL;
    srand(25);

    if(num>0)
    {
        double area=compute_surface_area(m);
        dense=sqrt(area/(double)num);
    }

    for(auto cit=m->c_begin(m);cit!=m->c_end(m);cit++)
    {
        sampling_from_one_cell(m,quote(cit),&re,dense );
    }

    return re;
}


static double** get_half_section_polygonxyz(double** poly,int len,int *len1 ,
    double axisz,int gq,int dim,int* tuihua)
{
    *len1=0;*tuihua=0;
    if(poly==NULL||len<=0)
    {
        return NULL;
    }

    int max_inter=10;

    double p[3]={0},pp[3]={0};

    int * ranse1=(int*)malloc(sizeof(int)*(len+max_inter));
    memset(ranse1,0,sizeof(int)*(len+max_inter));

    int ii=-1;
    double ** re=NULL;

    if(gq)
    {
        for(int i=0;i<len;i++)
        {
        //printf("temp_poly 1:%lf\n",temp_poly[i][1]);
            if(poly[i][dim]>axisz+1e-6)
            {
                ranse1[i]=1;
                ii=i;
            }
            else if(poly[i][dim]< axisz-1e-6)
            {
                ranse1[i]=-1;
            }
            else
            {
                *tuihua=*tuihua+1;
            }
        }
    }
    else
    {
        for(int i=0;i<len;i++)
        {
        //printf("temp_poly 1:%lf\n",temp_poly[i][1]);
            if(poly[i][dim]<axisz+1e-6)
            {
                ranse1[i]=1;
                ii=i;
            }
            else if(poly[i][dim]>axisz-1e-6)
            {
                ranse1[i]=-1;
            }
            else
            {
                *tuihua=*tuihua+1;
            }
        }

    }


    if(ii!=-1)
    {
        *len1=0;
        re=(double**)malloc(sizeof(double*)*(len+max_inter));
        memset(re,0,sizeof(double*)*(len+max_inter));
        int flag=1;
        for(int i=0;i<len;i++)
        {
            // printf("ranse:%d\n",ranse1[(i+ii)%len]);
            if(ranse1[(i+ii)%len]>0||(ranse1[(i+ii)%len]==0&&flag==1) )
            {
                re[*len1] =(double*)malloc(sizeof(double)*3);
                re[*len1][0]=poly[(i+ii)%len][0];re[*len1][1]=poly[(i+ii)%len ][1];re[*len1][2]=poly[(i+ii)%len][2];
                *len1=*len1+1;
            }
            if(ranse1[(i+ii)%len]>0){flag=1;}else if(ranse1[(i+ii)%len]<0){flag=-1;}

            if(( ranse1[(i+ii)%len]>0&& ranse1[(i+ii+1)%len]<0)||
                ( ranse1[(i+ii)%len]<0&& ranse1[(i+ii+1)%len]>0) )
            {

                double lambda=(poly[(i+ii+1)%len][dim]-axisz)/(poly[(i+ii+1)%len][dim]-poly[(i+ii)%len][dim]);
                p[0]=lambda*poly[(i+ii)%len][0]+(1-lambda)*poly[(i+ii+1)%len][0];
                p[1]=lambda*poly[(i+ii)%len][1]+(1-lambda)*poly[(i+ii+1)%len][1];
                p[2]=lambda*poly[(i+ii)%len][2]+(1-lambda)*poly[(i+ii+1)%len][2];

                re[*len1] =(double*)malloc(sizeof(double)*3);
                re[*len1][0]=p[0];re[*len1][1]=p[1];re[*len1][2]=p[2];
                *len1=*len1+1;
            }
        }
    }


    SAFE_FREE(ranse1);

    return re;
}

static int is_polygon_intersecting_with_bounding_box(double* min,double* max,double** poly,int len)
{
    int tuihua=0;
    int len1=len;

    double** re=(double**)malloc(sizeof(double*)*len);
    for(int i=0;i<len;i++)
    {
        re[i]=(double*)malloc(sizeof(double)*3);
        re[i][0]=poly[i][0];re[i][1]=poly[i][1];re[i][2]=poly[i][2];
    }
    for(int i=0;i<3;i++)
    {
        int len2=0;
        double** poly1= get_half_section_polygonxyz(re,len1,&len2,
            min[i],1,i,&tuihua);
        for(int i=0;i<len1;i++)
        {
            free(re[i]);
        }
        free(re);
        re=poly1;len1=len2;
        if(re==NULL){return 0;}
    }
    for(int i=0;i<3;i++)
    {
        int len2=0;
        double** poly1= get_half_section_polygonxyz(re,len1,&len2,
            max[i],0,i,&tuihua);
        for(int i=0;i<len1;i++)
        {
            free(re[i]);
        }
        free(re);
        re=poly1;len1=len2;
        if(re==NULL){return 0;}
    }

    for(int i=0;i<len1;i++)
    {
        free(re[i]);
    }
    free(re);
    return 1;
}

int oc_node_divide_one_leaf_on_cells(OC_Node* ocn)
{
    if(ocn==NULL)
    {return 0;}
    if(ocn->children[0]!=NULL)
    {
        printf("oc_node_divide_one_leaf is not leaf\n");
        return 0;
    }
    double steps[3]={0,0,0};
    for(int i=0;i<3;i++)
    {
        steps[i]=(ocn->loc_max[i]-ocn->loc_min[i])/2.0;
    }

    if(steps[0]<=0||steps[1]<=0||steps[2]<=0)
    {
        printf("steps is zero\n");
        return 0;
    }

    for(int i=0;i<2;i++)
    {
        for(int j=0;j<2;j++)
        {
            for(int k=0;k<2;k++)
            {
                ocn->children[i*4+j*2+k]=(OC_Node*)malloc(sizeof(OC_Node));
                oc_node_init(ocn->children[i*4+j*2+k]);
                ocn->children[i*4+j*2+k]->depth=ocn->depth+1;
                ocn->children[i*4+j*2+k]->parent=ocn;
                ocn->children[i*4+2*j+k]->loc_min[0]=ocn->loc_min[0]+steps[0]*i;
                ocn->children[i*4+2*j+k]->loc_min[1]=ocn->loc_min[1]+steps[1]*j;
                ocn->children[i*4+2*j+k]->loc_min[2]=ocn->loc_min[2]+steps[2]*k;
                ocn->children[i*4+2*j+k]->loc_max[0]=ocn->children[i*4+2*j+k]->loc_min[0]+steps[0];
                ocn->children[i*4+2*j+k]->loc_max[1]=ocn->children[i*4+2*j+k]->loc_min[1]+steps[1];
                ocn->children[i*4+2*j+k]->loc_max[2]=ocn->children[i*4+2*j+k]->loc_min[2]+steps[2];
            }
        }
    }

    // for(Node* nit=(Node*)(ocn->value);nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     template_v* vv=(template_v*)(nit->value);
    //     int i=0,j=0,k=0;
    //     i=(int)((vv->point[0]-ocn->loc_min[0])/steps[0]);
    //     j=(int)((vv->point[1]-ocn->loc_min[1])/steps[1]);
    //     k=(int)((vv->point[2]-ocn->loc_min[2])/steps[2]);
    //     i>=2?(i=1):(i<0?(i=0):(i=i));
    //     j>=2?(j=1):(j<0?(j=0):(j=j));
    //     k>=2?(k=1):(k<0?(k=0):(k=k));
    //     Node* n=(Node*)(ocn->children[i*4+j*2+k]->value);
    //     n=node_overlying(n,vv);
    //     ocn->children[i*4+2*j+k]->value=n;
    // }
    // free_node((Node*)(ocn->value));
    // ocn->value=NULL;
    int re=0;
    // for(int i=0;i<8;i++)
    // {   if(ocn->children[i]->value!=NULL)
    //     {re++;}
    // }
    return re;
}


static int  is_legal_swap_face(template_f* f)
{
    if(f->halffaces[0]->cell==NULL||f->halffaces[1]->cell==NULL)
    {return 0;}
    template_c* c1=f->halffaces[0]->cell;
    template_c* c2=f->halffaces[1]->cell;
    template_v*v1=NULL,*v2=NULL;
    for(int i=0;i<c1->vertices_size;i++)
    {
        int flag=1;
        for(int j=0;j<f->vertices_size;j++)
        {
            if(f->vertices[j]==c1->vertices[i])
            {
                flag=0;
                break;
            }
        }
        if(flag==1)
        {
            v1=c1->vertices[i];
            break;
        }
    }
    for(int i=0;i<c2->vertices_size;i++)
    {
        int flag=1;
        for(int j=0;j<f->vertices_size;j++)
        {
            if(f->vertices[j]==c2->vertices[i])
            {
                flag=0;
                break;
            }
        }
        if(flag==1)
        {
            v2=c2->vertices[i];
            break;
        }
    }
    // for(Node* nit=v1->cells;nit!=NULL;nit=(Node*)(nit->Next))
    for(int i=0;i<v1->cellsn.size;i++)
    {
        // for(Node* nit1=v2->cells;nit1!=NULL;nit1=(Node*)(nit1->Next))
        for(int j=0;j<v2->cellsn.size;j++)
        {
            if(v1->cellsn.extra_data[i]==v2->cellsn.extra_data[j])
            {
                return 0;
            }
        }
    }
    return 1;

}
Node* swap_face(Mesh* m,template_f* f)
{
    if(f==NULL||f->halffaces[0]->cell==NULL||f->halffaces[1]->cell==NULL)
    {return NULL;}
    template_c* c1=f->halffaces[0]->cell;
    template_c* c2=f->halffaces[1]->cell;
    template_v*v1=NULL,*v2=NULL;
    for(int i=0;i<c1->vertices_size;i++)
    {
        int flag=1;
        for(int j=0;j<f->vertices_size;j++)
        {
            if(f->vertices[j]==c1->vertices[i])
            {
                flag=0;
                break;
            }
        }
        if(flag==1)
        {
            v1=c1->vertices[i];
            break;
        }
    }
    for(int i=0;i<c2->vertices_size;i++)
    {
        int flag=1;
        for(int j=0;j<f->vertices_size;j++)
        {
            if(f->vertices[j]==c2->vertices[i])
            {
                flag=0;
                break;
            }
        }
        if(flag==1)
        {
            v2=c2->vertices[i];
            break;
        }
    }

    // for(Node* nit=v1->cells;nit!=NULL;nit=(Node*)(nit->Next))
    for(auto vcit=m->vc_begin(m,*v1);vcit!=m->vc_end(m,*v1);vcit++)
    {
        // for(Node* nit1=v2->cells;nit1!=NULL;nit1=(Node*)(nit1->Next))
        for(auto vcit1=m->vc_begin(m,*v2);vcit1!=m->vc_end(m,*v2);vcit1++)
        {
            if(quote(vcit)==quote(vcit1))
            {
                return NULL;
            }
        }
    }
    int size=c2->vertices_size;
    Node* nodec=NULL;
    for(auto chfit= m->chf_begin(m,*c2);chfit!=m->chf_end(m,*c2);chfit++)
    {
        if(quote(chfit)->face==f){continue;}
        template_v** vs=(template_v**)malloc(sizeof(template_v*)*size);
        memmove(vs,quote(chfit)->vertices,sizeof(template_v*)*quote(chfit)->vertices_size);
        vs[quote(chfit)->vertices_size]=v1;
        nodec=node_overlying(nodec,vs);

    }
    m->delete_cell(m,*c1,true);m->delete_cell(m,*c2,true);
    Node* re=NULL;
    for(Node* nit=nodec;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v**vs=(template_v**)(nit->value);
        template_c*nc=m->create_cellv(m,vs,size);
        re=node_overlying(re,nc);
        free(vs);
    }
    //free_node_value(nodec);
    free_node(nodec);
    return re;
}


static inline int is_edge_legacy_to_collapse(Mesh* m,template_f* f)
{
    LB_Vector_voidptr vec=Mesh_vv_beginn(m, *(f->vertices[0]));
    // LB_Vector_voidptr* vecp1=&vec;
    for(int i=1;i<f->vertices_size;i++)
    {
        if(vec.size==0){break;}
        // vec.intersection();
        LB_Vector_voidptr vec2=Mesh_vv_beginn(m, *(f->vertices[i]));
        LB_Vector_voidptr vec1= vec.intersection(&vec,&vec2);
        vec.clear(&vec);
        vec2.clear(&vec2);
        vec=vec1;

    }
    int size= vec.size;
    // free_node(n1);
    vec.clear(&vec);
    return size;
}

static int face_is_inter_with_face(template_f* f,template_hf* hf)
{
    if(hf->face==f)
    {
        return 1;
    }
    for(int i=0;i<f->vertices_size;i++)
    {
        for(int j=0;j<hf->vertices_size;j++)
        {
            if(f->vertices[i]==hf->vertices[j])
            {
                return 1;
            }
        }
    }
    return 0;
}
static void quchu_filter_halfface(Int_RB_Tree* tree_hf,template_f*f)
{
    Node* n=NULL;
    for(auto it=tree_hf->begin(tree_hf);it.it!=NULL;it++)
    {
        template_hf* hf=(template_hf*)(it.second);
        if(face_is_inter_with_face(f,hf)==1)
        {
            n=node_overlying(n,hf);
        }
    }
    for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf* hf=(template_hf*)(nit->value);
        tree_hf->erase(tree_hf,hf->id);

    }
    free_node(n);
}


template_v* collapse_edge(Mesh* m,template_f* f ,double* p)
{
    if(is_edge_legacy_to_collapse(m,f)>2)
    {
        return NULL;
    }

    // if(m->vertex_is_boundary(m,*(f->vertices[0]))||m->vertex_is_boundary(m,*(f->vertices[1])))
    // {
    //    // printf("cuowu1\n");
    //     return NULL;
    // }
    Int_RB_Tree* tree_c=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_c);
    for(int i=0;i<f->vertices_size;i++)
    {
        // for(Node* nit=f->vertices[i]->cells;nit!=NULL;nit=(Node*)(nit->Next))
        template_v*v=f->vertices[i];
        for(auto vcit=m->vc_begin(m,*v);vcit!=m->vc_end(m,*v);vcit++)
        {
            template_c* c=quote(vcit);
            // template_c* c=(template_c*)(nit->value);
            tree_c->insert(tree_c,c->id,c);
        }
    }
    //Node* bound_hfs=is_edge_legacy_to_collapse(m,f,tree_c);
    Int_RB_Tree* tree_hf=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_hf);
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
        {
            template_hf* hf=quote(chfit);

            if(tree_hf->find(tree_hf,hf->face->id)==NULL )
            {tree_hf->insert(tree_hf,hf->face->id,hf);}
            else
            {tree_hf->erase(tree_hf,hf->face->id);}
        }
    }

    for(int i=0;i<f->vertices_size;i++)
    {
        template_v*v=f->vertices[i];
        for(auto vfit=m->vf_begin(m,*v);vfit!=m->vf_end(m,*v);vfit++)
        // for(Node* nit=f->vertices[i]->faces;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_f* f=quote(vfit);
            // template_f* f=(template_f*)(nit->value);
            tree_hf->erase(tree_hf,f->id);
        }
    }

    int cols=f->vertices[0]->point_size;
    int rows=f->vertices_size+1;
    template_v*re=m->create_vertexv(m,p,cols);
    Node* node_hf=NULL;
    for(auto it=tree_hf->begin(tree_hf);it.it!=NULL;it++)
    {
        template_hf* hf=(template_hf*)(it.second);
        template_v** vs_value=(template_v**)malloc(sizeof(template_v*)*(hf->vertices_size+1) );
        memmove(vs_value,hf->vertices,sizeof(template_v*)*hf->vertices_size);
        vs_value[hf->vertices_size]=re;
        node_hf=node_overlying(node_hf,vs_value);
        //it.it->data=vs_value;
    }
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        m->delete_cell(m,*c,true);
    }
    for(Node* nit=node_hf;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v**vs_value=(template_v**)(nit->value);
        m->create_cellv(m,vs_value,rows);
        free(vs_value);
    }
    free_node(node_hf);
    //int_rb_tree_free_value(tree_hf);
    int_rb_tree_free(tree_hf);
    int_rb_tree_free(tree_c);
    return re;

}
static int compute_cell_normal(template_c*c,double*nor)
{
    if(c==NULL)
    {
        return 0;
    }

    template_hf* hf=(template_hf*)(c->halffaces[0]);
    template_v* v=NULL;
    for(int i=0;i<c->vertices_size;i++)
    {
        if(c->vertices[i]!=hf->vertices[0]&&c->vertices[i]!=hf->vertices[1])
        {
            v=c->vertices[i];
            break;
        }
    }

    if(v==NULL)
    {return 0;}

    double p1[3]={hf->vertices[0]->point[0]-v->point[0],hf->vertices[0]->point[1]-v->point[1],
        hf->vertices[0]->point[2]-v->point[2]};

    double p2[3]={hf->vertices[1]->point[0]-v->point[0],hf->vertices[1]->point[1]-v->point[1],
        hf->vertices[1]->point[2]-v->point[2]};

    out_productn(p1,p2,nor);
    normalize(nor,3);
    return 1;
    // double*re=out_product(p1,p2);
    // normalize(re,3);
    // return re;
}

static inline int compute_triangle_normal(double* q1,
    double* q2,double* q3,double*nor)
{
    double p1[3]={q2[0]-q1[0],q2[1]-q1[1],
        q2[2]-q1[2]};
    double p2[3]={q3[0]-q1[0],q3[1]-q1[1],
        q3[2]-q1[2]};
    out_productn(p1,p2,nor);
    normalize(nor,3);
    return 1;
}

static int is_face_fold(template_f*f,double threshold)
{
    if(f->halffaces[0]->cell==NULL||f->halffaces[1]->cell==NULL)
    {
        return 0;
    }
    double nor1[3]={0},nor2[3]={0};
    compute_cell_normal( f->halffaces[0]->cell,nor1);
    compute_cell_normal( f->halffaces[1]->cell,nor2);
    if(inner_product(nor1, nor2,3)<threshold)
    {
        return 1;
    }
    return 0;
}
static int is_face_swap(template_f*f,double threshold)
{
    template_c* c1=f->halffaces[0]->cell;
    template_c* c2=f->halffaces[1]->cell;
    if(c1==NULL||c2==NULL){return 0;}
    template_v*vs[4]={NULL};
    vs[0]=f->halffaces[0]->vertices[0]; vs[2]= f->halffaces[0]->vertices[1];
    for(int i=0;i<3;i++)
    {
        template_v* v=c1->vertices[i];
        if(v!=vs[0]&&v!=vs[2])
        {
            vs[3]=v;
            break;
        }
    }
    for(int i=0;i<3;i++)
    {
        template_v* v=c2->vertices[i];
        if(v!=vs[0]&&v!=vs[2])
        {
            vs[1]=v;
            break;
        }
    }
    double nor1[3]={0},nor2[3]={0};
    template_v* vs1[3]={vs[3],vs[0],vs[1]};
    compute_triangle_normal(vs1[0]->point,vs1[1]->point,vs1[2]->point,
        nor1);

    vs1[0]=vs[1];vs1[1]=vs[2];vs1[2]=vs[3];
    compute_triangle_normal(vs1[0]->point,vs1[1]->point,vs1[2]->point,
        nor2);
    if(inner_product(nor1, nor2, 3)<threshold)
    {
        return 0;
    }

    return 1;
}

void  adjust_fold_cells_swap_face(Mesh*m,double threshold)
{
    // Node* node_faces=NULL;
    Int_RB_Tree* treef=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(treef);
    for(auto fit=m->f_begin(m);fit!=m->f_end(m);fit++)
    {
        if(is_face_fold(quote(fit),threshold)==1)
        {
            treef->insert(treef,quote(fit)->id,quote(fit));
            // node_faces=node_overlying(node_faces, quote(fit));
        }
    }

    // for(Node* nit=node_faces;nit!=NULL;nit=(Node*)(nit->Next))
    int num=0;
    while(num<10)
    {
        printf("once fold cells：%d %d\n",num,treef->size);
        Int_RB_Tree*temp_tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
        int_rb_tree_init(temp_tree);

        for(auto it=treef->begin(treef);it.it!=NULL;it++)
        {
            template_f*f=m->get_facep(m,it.first);
            if(f==NULL){
                continue;}
            if(is_face_fold(f,threshold)==0)
            {
                continue;
            }
            temp_tree->insert(temp_tree,f->id,f);
            if(is_face_swap(f,threshold)==0)
            {
                continue;
            }
            Node* cellsn=swap_face(m, f);
            for(Node* nit=cellsn;nit!=NULL;nit=(Node*)(nit->Next))
            {
                template_c*c=(template_c*)(nit->value);
                for(int i=0;i<c->halffaces_size;i++)
                {
                    temp_tree->insert(temp_tree,
                        c->halffaces[i]->face->id,
                        c->halffaces[i]->face);
                }
            }
            free_node(cellsn);
        }
        int_rb_tree_free(treef);
        treef=temp_tree;
        num++;
    }

    int_rb_tree_free(treef);
    // free_node(node_faces);
}

int is_cells_simple_connected_on_surface(Mesh* m,Int_RB_Tree* cells)
{
    Int_RB_Tree* treev=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(treev);

    for(auto it=cells->begin(cells);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        for(int i=0;i<c->vertices_size;i++)
        {
            treev->insert(treev,c->vertices[i]->id,c->vertices[i]);
        }
    }

    Int_RB_Tree* treev1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(treev1);

    for(auto it=treev->begin(treev);it.it!=NULL;it++)
    {
        template_v* v=(template_v*)(it.second);
        if(m->vertex_is_boundary(m,*v))
        {
            continue;
        }
        for(int i=0;i<v->cellsn.size;i++)
        {
            template_c*c=(template_c*)(
                *lb_vector1_at_voidptr(&(v->cellsn), i));
            if(cells->find(cells,c->id)==NULL)
            {
                treev1->insert(treev1,v->id,v);
                break;
            }
        }
    }

    LB_Vector_voidptr vec,vec1;
    lb_vector_init_voidptr(&vec);
    lb_vector_init_voidptr(&vec1);

    template_v* v=(template_v*)(treev1->begin(treev1).second);
    treev1->erase(treev1,v->id);
    vec.push_back(&vec,v);
    LB_Vector_voidptr *vecp=&vec,*vecp1=&vec1;
    int flag=1;
    while(flag==1)
    {
        flag=0;
        for(int i=0;i<vecp->size;i++)
        {
            template_v*vv=(template_v*)(*(vecp->at(vecp,i)));
            for(int j=0;j<vv->cellsn.size;j++)
            {
                template_c*c=(template_c*)(*lb_vector1_at_voidptr(&(vv->cellsn), j));
                if(cells->find(cells,c->id)!=NULL)
                {continue;}
                for(int k=0;k<c->vertices_size;k++)
                {
                    template_v*vvv=c->vertices[k];
                    if(treev1->find(treev1,vvv->id)!=NULL)
                    {
                        treev1->erase(treev1,vvv->id);
                        vecp1->push_back(vecp1,vvv);
                        flag=1;
                    }
                }

            }
            // LB_Vector_voidptr vecc= Mesh_vv_beginn(m,*vv);
            // for(int j=0;j<vecc.size;j++)
            // {
            //     template_v*vvv=(template_v*)(*(vecc.at(&vecc,j)));
            //     if(treev1->find(treev1,vvv->id)!=NULL)
            //     {
            //         treev1->erase(treev1,vvv->id);
            //         vecp1->push_back(vecp1,vvv);
            //         flag=1;
            //     }
            // }
            // vecc.clear(&vecc);
        }
        vecp->clear(vecp);
        LB_Vector_voidptr*temp_vec=vecp;
        vecp=vecp1;
        vecp1=temp_vec;
    }
    int re=0;
    if(treev1->size<=0){re=1;}

    vec.clear(&vec);

    vec1.clear(&vec1);
    int_rb_tree_free(treev1);

    int_rb_tree_free(treev);

    return re;
}
// template_v* collapse_edgen(Mesh* m,template_f* f ,double* p,double* normal)
// {
//     if(is_edge_legacy_to_collapse(m,f)>2)
//     {
//         return NULL;
//     }

//     // if(m->vertex_is_boundary(m,*(f->vertices[0]))||m->vertex_is_boundary(m,*(f->vertices[1])))
//     // {
//     //    // printf("cuowu1\n");
//     //     return NULL;
//     // }
//     Int_RB_Tree* tree_c=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
//     int_rb_tree_init(tree_c);
//     for(int i=0;i<f->vertices_size;i++)
//     {
//         for(Node* nit=f->vertices[i]->cells;nit!=NULL;nit=(Node*)(nit->Next))
//         {
//             template_c* c=(template_c*)(nit->value);
//             tree_c->insert(tree_c,c->id,c);
//         }
//     }
//     //Node* bound_hfs=is_edge_legacy_to_collapse(m,f,tree_c);
//     Int_RB_Tree* tree_hf=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
//     int_rb_tree_init(tree_hf);
//     for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
//     {
//         template_c* c=(template_c*)(it.second);
//         for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
//         {
//             template_hf* hf=quote(chfit);

//             if(tree_hf->find(tree_hf,hf->face->id)==NULL )
//             {tree_hf->insert(tree_hf,hf->face->id,hf);}
//             else
//             {tree_hf->erase(tree_hf,hf->face->id);}
//         }
//     }

//     for(int i=0;i<f->vertices_size;i++)
//     {
//         for(Node* nit=f->vertices[i]->faces;nit!=NULL;nit=(Node*)(nit->Next))
//         {
//             template_f* f=(template_f*)(nit->value);
//             tree_hf->erase(tree_hf,f->id);
//         }
//     }

//     int cols=f->vertices[0]->point_size;
//     int rows=f->vertices_size+1;
//     template_v*re=m->create_vertexv(m,p,cols);
//     Node* node_hf=NULL;
//     for(auto it=tree_hf->begin(tree_hf);it.it!=NULL;it++)
//     {
//         template_hf* hf=(template_hf*)(it.second);
//         template_v** vs_value=(template_v**)malloc(sizeof(template_v*)*(hf->vertices_size+1) );
//         memmove(vs_value,hf->vertices,sizeof(template_v*)*hf->vertices_size);
//         vs_value[hf->vertices_size]=re;
//         node_hf=node_overlying(node_hf,vs_value);
//         //it.it->data=vs_value;
//     }
//     for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
//     {
//         template_c*c=(template_c*)(it.second);
//         m->delete_cell(m,*c,true);
//     }
//     for(Node* nit=node_hf;nit!=NULL;nit=(Node*)(nit->Next))
//     {
//         template_v**vs_value=(template_v**)(nit->value);
//         m->create_cellv(m,vs_value,rows);
//         free(vs_value);
//     }
//     free_node(node_hf);
//     //int_rb_tree_free_value(tree_hf);
//     int_rb_tree_free(tree_hf);
//     int_rb_tree_free(tree_c);
//     return re;
// }

// static int libcell_tools_two_halfaces_connected(template_hf* hf,template_hf* hf1)
// {
//     for(int i=0;i<hf->vertices_size;i++)
//     {
//         for(int j=hf1->vertices_size-1;j>=0;j--)
//         {
//             if(hf->vertices[i]==hf1->vertices[j])
//             {
//                 return 1;
//             }
//         }
//     }
//     return 0;
// }

// static Node*  is_edge_legacy_to_collapse(Mesh* m,template_f* f,Int_RB_Tree* tree_c)
// {

//     Int_RB_Tree* tree_hf=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
//     int_rb_tree_init(tree_hf);
//     for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
//     {
//         template_c* c=(template_c*)(it.second);
//         for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++)
//         {
//             template_hf* hf=quote(chfit);
//             if(tree_hf->find(tree_hf,hf->id)==NULL)
//             {
//                 tree_hf->insert(tree_hf,hf->id,hf);
//             }
//             else
//             {
//                 tree_hf->erase(tree_hf,hf->id);
//             }
//         }
//     }
//     Node* re=NULL;
//     if(tree_hf->size>0)
//     {
//         template_hf* hf=(template_hf*)((tree_hf->begin(tree_hf)).second);
//         re=node_overlying(re,hf);
//         tree_hf->erase(tree_hf,hf->id);
//     }

//     int flag=1;
//     while(flag)
//     {
//         flag=0;
//         template_hf* hf=(template_hf*)(re->value);
//         for(int i=0;i<hf->vertices_size;i++)
//         {
//             for(auto it= tree_hf->begin(tree_hf);it.it!=NULL;it++)
//             {
//                 template_hf* hf1=(template_hf*)(it.second);
//                 for(int j=hf1->vertices_size-1;j>=0;j--)
//                 {
//                     if(hf1->vertices[j]==hf->vertices[i])
//                     {
//                         flag=1;
//                         break;
//                     }
//                 }
//                 if(flag==1)
//                 {
//                     flag=1;
//                     re=node_overlying(re,hf1);
//                     tree_hf->erase(tree_hf,hf1->id);
//                     break;
//                 }
//             }
//             if(flag){break;}
//         }
//     }
//     if(tree_hf->size>0)
//     {
//         free_node(re);
//         re=NULL;
//     }
//     int_rb_tree_free(tree_hf);

//     return re;
// }
